Otimização de rotas:

Utilizando OR-Tools e OpenStreetMap

Thiago Pires | IBM

Thiago Pires

Formação:

  • Estatística (UERJ)
  • MSc. Epidemiologia (ENSP/FIOCRUZ)
  • DSc. Engenharia Biomédica (COPPE/UFRJ)


Experiência profissional:

  • Pesquisador (FIOCRUZ)
  • Professor (UFF)
  • Estatístico (FGV)
  • Cientista de Dados (IBM)

Problema do caxeiro viajante

O problema do caixeiro-viajante envolve encontrar a rota mais curta para visitar várias cidades uma vez cada, retornando à cidade de partida. É inspirado na otimização das entregas de vendedores para economizar tempo, custos de transporte e combustível.

Problema de roteamento do veículo

No Problema de roteamento de veículo (VRP, na sigla em inglês), o objetivo é encontrar as melhores rotas para vários veículos que visitam um conjunto de locais. Quando há apenas um veículo, ele é reduzido para o mesmo problema do caxeiro viajante.


Restrições

  • Restrições de capacidade: os veículos precisam pegar os itens em cada local que visitarem, mas têm uma capacidade máxima de transporte.

  • Janelas de tempo: cada local precisa ser visitado em uma janela de tempo específica.

Problema de roteamento de veículo com janelas de tempo

O VRPTW (sigla em inglês) é uma variação do VRP clássico, e pode ser descrito como um problema no qual uma frota de veículos que tem como origem um depósito central deve atender um conjunto de clientes que possuem suas próprias demandas.


Cada cliente deve ser atendido dentro de sua janela de tempo, que são horários predefinidos pelos mesmos para receber suas cargas.

Exemplo de otimização

Variáveis do problema

  • 3 Técnicos para atender uma região
  • 10 chamados para serem atendidos
  • As localizações do chamado tenho informação de CEP
  • Cada chamado tem um Service Level Agreement de atendimento. Assim, em determinada janela de tempo precisa-se iniciar o atendimento

Dados do problema

import pandas as pd

chamado = [str(i).rjust(2, '0') for i in range(1, 11)]
cep = ['22040010', '22010030', '22231200', '20511270', '20510110',
       '20560010', '20715040', '22250906', '20250430', '20261005']

chamados = pd.DataFrame({'chamado': chamado, 'cep': cep})
chamados.head(10)
chamado cep
0 01 22040010
1 02 22010030
2 03 22231200
3 04 20511270
4 05 20510110
5 06 20560010
6 07 20715040
7 08 22250906
8 09 20250430
9 10 20261005

Geolocalização dos chamados

Na otimização iremos precisar calcular as distâncias entre todos os chamados para calcular uma matriz de custo. Então precisamos calcular a geolocalizações.


import requests
from geopy.geocoders import Nominatim

geolocator = Nominatim(user_agent="th1460")

def cep2coo(cep):
    response = requests.get(f'https://viacep.com.br/ws/{cep}/json/')
    results = response.json()
    try:
        query = f"{results['logradouro']}, {results['localidade']} {results['uf']}"
        data = geolocator.geocode(query)
        return data.point.longitude, data.point.latitude
    except:
        return ()

chamados['coo'] = chamados['cep'].apply(cep2coo)

Persistindo os dados

import duckdb

con = duckdb.connect('../data/chamados.duckdb')
con.sql('CREATE TABLE chamados AS SELECT * FROM chamados')
con.close()

Lendo os dados

import duckdb

con = duckdb.connect('../data/chamados.duckdb')
chamados = con.sql('SELECT * FROM chamados;').fetchdf()
con.close()

chamados.head(5)
chamado cep coo
0 01 22040010 [-43.1835746, -22.9672354]
1 02 22010030 [-43.1667647, -22.9622599]
2 03 22231200 [-43.1851476, -22.9355907]
3 04 20511270 [-43.2412904, -22.9286744]
4 05 20510110 [-43.2416345, -22.9258574]

Definir limites (bouding box) da região

Esta definição é importante para limitar a região a ser analisada.

import osmnx as ox

lng = list(zip(*chamados.coo))[0]
lat = list(zip(*chamados.coo))[1]
place = [min(lat) + .1, max(lat) - .1, min(lng) - .1, max(lng) + .1]
place
[-22.8672354, -23.0112068, -43.364659, -43.0667647]

Região no mapa

import folium
m = folium.Map(
    location=[-22.939111, -43.219046],
    zoom_start=11
)

for site in chamados.iterrows():
    folium.Marker(list(reversed(site[1]['coo'])), 
    popup=site[1]['chamado']).add_to(m)

bounds = [(min(lat) - .05, min(lng) - .05), 
          (max(lat) + .05, max(lng) + .05)]
folium.Rectangle(bounds=bounds, color='#526D82').add_to(m)
m
Make this Notebook Trusted to load map: File -> Trust Notebook

Contruir grafo da região

graph = ox.graph_from_bbox(*place, simplify=True, network_type='drive')
nodes, streets = ox.graph_to_gdfs(graph)

Nodes

nodes.head()
y x street_count highway ref geometry
osmid
25038483 -22.980949 -43.206859 3 NaN NaN POINT (-43.20686 -22.98095)
25038484 -22.981045 -43.204756 5 NaN NaN POINT (-43.20476 -22.98105)
25038485 -22.981160 -43.202639 4 NaN NaN POINT (-43.20264 -22.98116)
25038487 -22.981267 -43.200507 4 NaN NaN POINT (-43.20051 -22.98127)
25038488 -22.981470 -43.206925 4 NaN NaN POINT (-43.20693 -22.98147)

Streets

streets.head(3)
osmid oneway lanes name highway reversed length geometry maxspeed bridge ref tunnel junction access width area
u v key
25038483 25038488 0 4217297 True 2 Rua Maria Quitéria residential False 58.357 LINESTRING (-43.20686 -22.98095, -43.20693 -22... NaN NaN NaN NaN NaN NaN NaN NaN
528897321 0 4217299 True 2 Rua Alberto de Campos residential False 47.768 LINESTRING (-43.20686 -22.98095, -43.20639 -22... NaN NaN NaN NaN NaN NaN NaN NaN
25038484 25038505 0 4217293 True 2 Rua Joana Angélica residential False 97.653 LINESTRING (-43.20476 -22.98105, -43.20472 -22... NaN NaN NaN NaN NaN NaN NaN NaN

Identificação dos nodes mais próximos de cada cep

tickets_osmid = ox.distance.nearest_nodes(graph, list(zip(*chamados.coo))[0], list(zip(*chamados.coo))[1])
chamados['osmid'] = tickets_osmid
chamados.head()
chamado cep coo osmid
0 01 22040010 [-43.1835746, -22.9672354] 148824293
1 02 22010030 [-43.1667647, -22.9622599] 153791875
2 03 22231200 [-43.1851476, -22.9355907] 331107980
3 04 20511270 [-43.2412904, -22.9286744] 1334341480
4 05 20510110 [-43.2416345, -22.9258574] 1752004535

Ligar nodes com os chamados

tickets_nodes = nodes[nodes.index.isin(tickets_osmid)]

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(10,10))
ax.set_facecolor('black')
streets.plot(ax=ax, color='gray', zorder=0, linewidth=1)
tickets_nodes.plot(ax=ax, color='red', zorder=5)
plt.grid()

Calcular matriz de custo

# impute speed on all edges
graph = ox.add_edge_speeds(graph)

# calculate travel time (seconds) for all edges
graph = ox.add_edge_travel_times(graph)

# update streets
streets = ox.graph_to_gdfs(graph, nodes=False)

Combinações de origem e destino

import numpy as np
import networkx as nx
from itertools import combinations

# route combinations
unique_tickets = np.unique(tickets_osmid)

route_combinations = list(combinations(unique_tickets, 2))

def cost_route(route):
    paths = nx.shortest_path(graph, route[0], route[1], weight='travel_time')
    return sum(streets.loc[paths,'travel_time'].to_list())/6

Calcular matriz de custo

total_travel_time = [cost_route(i) for i in route_combinations]
index = pd.MultiIndex.from_tuples(route_combinations, names=['orig', 'dest'])
df_total_travel_time = pd.Series(total_travel_time, index=index)
df_total_travel_time.drop_duplicates(inplace=True)
df_total_travel_time.head()
orig       dest     
148824293  149369298     87.566667
           153791875     43.083333
           331107980    126.216667
           506427643    224.916667
           837056712    303.466667
dtype: float64

Formatando a matriz

cost_matrix = np.zeros(shape=(len(unique_tickets), len(unique_tickets)))

for i, j in enumerate(unique_tickets):
    for k, l in enumerate(unique_tickets):

        try:
            cost_matrix[i, k] =  df_total_travel_time.loc[(j, l), ]
        except KeyError:
            cost_matrix[i, k] =  0

i_lower = np.tril_indices(len(unique_tickets), -1)
cost_matrix[i_lower] = cost_matrix.T[i_lower]
cost_matrix = cost_matrix.astype(int)
cost_matrix = np.insert(cost_matrix, 0, 0, axis=0)
cost_matrix = np.insert(cost_matrix, 0, 0, axis=1)
cost_matrix = cost_matrix.tolist()
cost_matrix
[[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
 [0, 0, 87, 43, 126, 224, 303, 200, 259, 268, 240],
 [0, 87, 0, 76, 44, 190, 256, 141, 215, 223, 194],
 [0, 43, 76, 0, 98, 245, 310, 195, 269, 277, 248],
 [0, 126, 44, 98, 0, 183, 248, 133, 207, 215, 186],
 [0, 224, 190, 245, 183, 0, 206, 128, 84, 93, 143],
 [0, 303, 256, 310, 248, 206, 0, 265, 177, 155, 96],
 [0, 200, 141, 195, 133, 128, 265, 0, 125, 133, 138],
 [0, 259, 215, 269, 207, 84, 177, 125, 0, 40, 110],
 [0, 268, 223, 277, 215, 93, 155, 133, 40, 0, 79],
 [0, 240, 194, 248, 186, 143, 96, 138, 110, 79, 0]]

Outros parâmetros

  • Pontos iniciais
starts = [6, 2, 3]
  • Pontos finais
ends = [0, 0, 0]
  • Janelas de tempo
time_windows = [(0, 0), (0, 300), (0, 420), (0, 500), (0, 340), (0, 220), (0, 332), (0, 234), (0, 423), (0, 350), (0, 630)]
  • Número de técnicos
num_engineers = 3

Otimização

Inputs

def create_data_model():
    """Stores the data for the problem."""
    data = {}
    data['time_matrix'] = cost_matrix
    data['time_windows'] = time_windows
    data['num_engineers'] = num_engineers
    data['starts'] = starts
    data['ends'] = ends
    return data

Otimização

from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp


def print_solution(data, manager, routing, solution):
    """Prints solution on console."""
    print(f'Objective: {solution.ObjectiveValue()}')
    time_dimension = routing.GetDimensionOrDie('Time')
    total_time = 0
    for engineer_id in range(data['num_engineers']):
        index = routing.Start(engineer_id)
        plan_output = 'Route for engineer {}:\n'.format(engineer_id)
        while not routing.IsEnd(index):
            time_var = time_dimension.CumulVar(index)
            plan_output += '{0} Time({1},{2}) -> '.format(
                manager.IndexToNode(index), solution.Min(time_var),
                solution.Max(time_var))
            index = solution.Value(routing.NextVar(index))
        time_var = time_dimension.CumulVar(index)
        plan_output += '{0} Time({1},{2})\n'.format(manager.IndexToNode(index),
                                                    solution.Min(time_var),
                                                    solution.Max(time_var))
        plan_output += 'Time of the route: {}min\n'.format(
            solution.Min(time_var))
        print(plan_output)
        total_time += solution.Min(time_var)
    print('Total time of all routes: {}min'.format(total_time))


def main():
    """Solve the VRP with time windows."""
    # Instantiate the data problem.
    data = create_data_model()

    # Create the routing index manager.
    manager = pywrapcp.RoutingIndexManager(len(data['time_matrix']),
                                           data['num_engineers'],
                                           data['starts'],
                                           data['ends'])

    # Create Routing Model.
    routing = pywrapcp.RoutingModel(manager)

    # Create and register a transit callback.
    def time_callback(from_index, to_index):
        """Returns the travel time between the two nodes."""
        # Convert from routing variable Index to time matrix NodeIndex.
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        return data['time_matrix'][from_node][to_node]

    transit_callback_index = routing.RegisterTransitCallback(time_callback)

    # Define cost of each arc.
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

    # Add Time Windows constraint.
    time = 'Time'
    routing.AddDimension(
        transit_callback_index,
        30,  # allow waiting time
        400,  # maximum time per engineer
        False,  # Don't force start cumul to zero.
        time)
    time_dimension = routing.GetDimensionOrDie(time)
    # Add time window constraints for each location except start locations.
    starts = data['starts'] + [0]
    for location_idx, time_window in enumerate(data['time_windows']):
        if location_idx in starts:
            continue
        index = manager.NodeToIndex(location_idx)
        time_dimension.CumulVar(index).SetRange(time_window[0], time_window[1])

    # Instantiate route start and end times to produce feasible times.
    for i in range(data['num_engineers']):
        routing.AddVariableMinimizedByFinalizer(
            time_dimension.CumulVar(routing.Start(i)))
        routing.AddVariableMinimizedByFinalizer(
            time_dimension.CumulVar(routing.End(i)))

    # Setting first solution heuristic.
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)

    # Solve the problem.
    solution = routing.SolveWithParameters(search_parameters)

    # Print solution on console.
    if solution:
        print_solution(data, manager, routing, solution)


if __name__ == '__main__':
    main()
Objective: 688
Route for engineer 0:
6 Time(0,0) -> 5 Time(206,206) -> 8 Time(290,290) -> 9 Time(330,330) -> 0 Time(330,330)
Time of the route: 330min

Route for engineer 1:
2 Time(0,0) -> 4 Time(44,44) -> 7 Time(177,177) -> 10 Time(315,315) -> 0 Time(315,315)
Time of the route: 315min

Route for engineer 2:
3 Time(0,0) -> 1 Time(43,43) -> 0 Time(43,43)
Time of the route: 43min

Total time of all routes: 688min

Otimização considerando skills

Inputs

def create_data_model():
    """Stores the data for the problem."""
    data = {}
    data['time_matrix'] = cost_matrix
    data['time_windows'] = time_windows
    data['skill'] = [
        [0, 1, 2, 3],
        [0, 1, 2, 3],
        [0, 1, 2, 3],
        [0, 1, 2, 3],
        [0, 1, 2, 3],
        [0, 1, 2, 3],
        [0, 1, 2, 3],
        [0, 1, 2, 3],
        [1, 2, 3],
        [0, 1, 2, 3]
    ]
    data['num_engineers'] = num_engineers
    data['starts'] = starts
    data['ends'] = ends
    return data

Otimização

from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp


def print_solution(data, manager, routing, solution):
    """Prints solution on console."""
    print(f'Objective: {solution.ObjectiveValue()}')
    time_dimension = routing.GetDimensionOrDie('Time')
    total_time = 0
    for engineer_id in range(data['num_engineers']):
        index = routing.Start(engineer_id)
        plan_output = 'Route for engineer {}:\n'.format(engineer_id)
        while not routing.IsEnd(index):
            time_var = time_dimension.CumulVar(index)
            plan_output += '{0} Time({1},{2}) -> '.format(
                manager.IndexToNode(index), solution.Min(time_var),
                solution.Max(time_var))
            index = solution.Value(routing.NextVar(index))
        time_var = time_dimension.CumulVar(index)
        plan_output += '{0} Time({1},{2})\n'.format(manager.IndexToNode(index),
                                                    solution.Min(time_var),
                                                    solution.Max(time_var))
        plan_output += 'Time of the route: {}min\n'.format(
            solution.Min(time_var))
        print(plan_output)
        total_time += solution.Min(time_var)
    print('Total time of all routes: {}min'.format(total_time))


def main():
    """Solve the VRP with time windows."""
    # Instantiate the data problem.
    data = create_data_model()

    # Create the routing index manager.
    manager = pywrapcp.RoutingIndexManager(len(data['time_matrix']),
                                           data['num_engineers'],
                                           data['starts'],
                                           data['ends'])

    # Create Routing Model.
    routing = pywrapcp.RoutingModel(manager)

    # Create and register a transit callback.
    def time_callback(from_index, to_index):
        """Returns the travel time between the two nodes."""
        # Convert from routing variable Index to time matrix NodeIndex.
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        return data['time_matrix'][from_node][to_node]

    transit_callback_index = routing.RegisterTransitCallback(time_callback)

    # Define cost of each arc.
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

    # Add Time Windows constraint.
    time = 'Time'
    routing.AddDimension(
        transit_callback_index,
        30,  # allow waiting time
        400,  # maximum time per engineer
        False,  # Don't force start cumul to zero.
        time)
    time_dimension = routing.GetDimensionOrDie(time)
    # Add time window constraints for each location except start locations.
    starts = data['starts'] + [0]
    for location_idx, time_window in enumerate(data['time_windows']):
        if location_idx in starts:
            continue
        index = manager.NodeToIndex(location_idx)
        time_dimension.CumulVar(index).SetRange(time_window[0], time_window[1])

    # Instantiate route start and end times to produce feasible times.
    for i in range(data['num_engineers']):
        routing.AddVariableMinimizedByFinalizer(
            time_dimension.CumulVar(routing.Start(i)))
        routing.AddVariableMinimizedByFinalizer(
            time_dimension.CumulVar(routing.End(i)))

    # Setting first solution heuristic.
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)

    # Constraint to skills (set engineers to each ticket)
    for ticket, engineers in enumerate(data['skill'], start=1):
        engineers.insert(0, -1)
        routing.VehicleVar(manager.NodeToIndex(ticket)).SetValues(engineers)

    # Solve the problem.
    solution = routing.SolveWithParameters(search_parameters)

    # Print solution on console.
    if solution:
        print_solution(data, manager, routing, solution)


if __name__ == '__main__':
    main()
Objective: 734
Route for engineer 0:
6 Time(0,0) -> 5 Time(206,206) -> 10 Time(349,349) -> 0 Time(349,349)
Time of the route: 349min

Route for engineer 1:
2 Time(0,0) -> 4 Time(44,44) -> 7 Time(177,177) -> 8 Time(302,302) -> 9 Time(342,342) -> 0 Time(342,342)
Time of the route: 342min

Route for engineer 2:
3 Time(0,0) -> 1 Time(43,43) -> 0 Time(43,43)
Time of the route: 43min

Total time of all routes: 734min